##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(ggpubr)
library(patchwork)
library("scales")

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"

    sample_list_file <- paste(work_dir,"/public_ref/combine/MutationInfo.combine.tsv",sep="")
    rsem_file <- "~/20220915_gastric_multiple/dna_combinePublic/mRNA/CombineTMM.DNAUse.NJMU_TCGA.MergeMutiSample.tsv"

    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/expression"
    gene <- "MUC6"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
out_path <- opt$out_path
rsem_file <- opt$rsem_file
gene <- opt$gene

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_expression <- data.frame(fread(rsem_file))

##########################################################################################
plotTpm <- function(tmp_dat = tmp_dat , out_name = out_name , title = title , from = from , y_max = y_tmm){

    if(from!="TCGA"){
        my_comparisons_1 <- list( 
            c(1, 2) , c(1, 3) , c(1, 4) , 
            c(2, 3) , c(2, 4) ,
            c(3, 4)  
        )
    }else{
        my_comparisons_1 <- list(  c(1, 2)  )
    }

    if(y_max > 2000){

        if(y_max > 10000){
            y_breaks <- 10000
        }else if(y_max > 5000){
            y_breaks <- 5000
        }else{
            y_breaks <- 2000
        }
        y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_max*2.5) , breaks = c( 0 , 500 , 1000 , 2000 , y_breaks ) , trans = sqrt_trans() )
    }else{
        y_lab_lim <- scale_y_continuous( limits = c(-0.1 ,y_max*1.5))
    }

    plot <- ggplot( tmp_dat , aes( x = Class , y = Tpm , color = Class ) ) +
        geom_boxplot(alpha =1 , outlier.color=NA , size = 0.9 , width = 0.6) +
        geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
        scale_fill_npg()+
        scale_color_npg()+
        xlab(NULL) +
        ylab("TMM")+
        theme_bw() +
        y_lab_lim +
        ggtitle(title) +
        stat_compare_means(comparisons = my_comparisons_1) +
        theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='none',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 12,color="black",face='bold',hjust=0.5,vjust=0.5),
            legend.text = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 10,color="black",face='bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 10,color="black",face='bold') ,
            axis.line = element_line(size = 0.5)) 

    return(plot)
}

describeGene <- function(dat_expression = dat_expression , image_path = image_path ){

    tmp_exp <- subset( dat_expression , gene_id == gene)

    sample <- sapply( strsplit(colnames(tmp_exp)[2:(ncol(tmp_exp))] , "_") , "[" , 1)
    class <- sapply( strsplit(colnames(tmp_exp)[2:(ncol(tmp_exp))] , "_") , "[" , 2)
    tpm <- as.numeric(tmp_exp[2:(ncol(tmp_exp))])
    tmp_dat <- data.frame( Sample = sample , Class = class , Tpm = tpm )
    tmp_dat$Class <- factor( tmp_dat$Class , levels = c( "Normal" , "IM" , "IGC" , "DGC" ) , order = T )

    y_tmm <- max(tmp_dat$Tpm) * 1.1

    from <- "All"
    tmp_dat_use <- tmp_dat
        title <- paste0(
        "Gene : " , gene , "\n" , from
    )
    out_name <- paste0( image_path , "/" , gene , ".normalize.oneImage.",from,".pdf" )
    plot <- plotTpm(tmp_dat = tmp_dat_use , out_name = out_name , title = title , from = from , y_max = y_tmm)
    ggsave(file=out_name,plot=plot,width=4,height=5)
    out_name <- paste0( image_path , "/" , gene , ".normalize.",from,".tsv" )
    write.table(tmp_dat_use , out_name , row.names = F , sep = "\t" , quote = F)

    ## NJMU和TCGA画在一张图上
    from <- "NJMU"
    index <- grep( "TCGA" , tmp_dat$Sample , invert = T)
    tmp_dat_use <- tmp_dat[index,]
    title <- paste0(
        "Gene : " , gene , "\n" , from
    )
    plot1 <- plotTpm(tmp_dat = tmp_dat_use , out_name = out_name , title = title , from = from , y_max = y_tmm)
    out_name <- paste0( image_path , "/" , gene , ".normalize.",from,".tsv" )
    write.table(tmp_dat_use , out_name , row.names = F , sep = "\t" , quote = F)

    from <- "TCGA"
    index <- grep( "TCGA" , tmp_dat$Sample )
    tmp_dat_use <- tmp_dat[index,]
        title <- paste0(
        "Gene : " , gene , "\n" , from
    )
    plot2 <- plotTpm(tmp_dat = tmp_dat_use , out_name = out_name , title = title , from = from , y_max = y_tmm)
    out_name <- paste0( image_path , "/" , gene , ".normalize.",from,".tsv" )
    write.table(tmp_dat_use , out_name , row.names = F , sep = "\t" , quote = F)

    out_name <- paste0( image_path , "/" , gene , ".normalize.oneImage.From.pdf" )
    ggsave(file=out_name,plot=plot1 + plot2,width=6,height=5)

}


##########################################################################################
## 分单个基因，看在一个人多个样本的改变情况
Hugo_Symbol <- gene
image_path <- out_path
dir.create(image_path , recursive = T)

describeGene(dat_expression = dat_expression , image_path = image_path )



